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The Fast Norm Vector Indicator (FN VI) method: A new dynamical 
parameter for detecting order and chaos in Hamiltonian systems 
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Abstract In the present article, we introduce and also de- 
ploy a new, simple, very fast and efficient method, the Fast 
Norm Vector Indicator (FNVI) in order to distinguish rapidly 
and with certainty between ordered and chaotic motion in 
Hamiltonian systems. This distinction is based on the dif- 
ferent behavior of the FNVI for the two cases: the indicator 
after a very short transient period of fluctuation displays a 
nearly constant value for regular orbits, while it continues 
to fluctuate significantly for chaotic orbits. In order to quan- 
tify the results obtained by the FNVI method, we establish 
the dFNVI, which is the quantified numerical version of the 
FNVI. A thorough study of the method's ability to achieve 
an early and clear detection of an orbit's behavior is pre- 
sented both in two and three degrees of freedom (2D and 
3D) Hamiltonians. Exploiting the advantages of the dFNVI 
method, we demonstrate how one can rapidly identify even 
tiny regions of order or chaos in the phase space of Hamil- 
tonian systems. The new method can also be applied in or- 
der to follow the time evolution of sticky orbits. A detailed 
comparison between the new FNVI method and some other 
well-known dynamical methods of chaos detection reveals 
the great efficiency and the leading role of this new dynam- 
ical indicator. 
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1 Introduction 

The issue of knowing whether the orbits of a dynamical 
system are ordered or chaotic is fundamental for the under- 
standing of the behavior of the system in an extended area 
of modem science. In the dissipative case, ffiis distinction 
is easily made as both types of motion are attracting. On 
ffie offier hand, in conservative systems, distinguishing be- 
tween regular and chaotic motion is often a very delicate and 
difficult issue (i.e. when the chaotic or ordered regions are 
small), especially in dynamical systems with many degrees 
of freedom, where one cannot easily visualize and interpret 
directly the behavior of the orbits. For this reason, it is of- 
ten of great importance to possess fast and accurate tools in 
order to determine if an orbit is ordered or chaotic, indepen- 
dent of the dimension of the phase space of the dynamical 
system. 

Let's recall and present some well-known, classical meth- 
ods that try to give an answer to the issue of determining the 
nature of an orbit. 

(a) The inspection of the consequents of an orbit on a 
Poincare surface of section (PSS). For 2D dynamical sys- 
tems this technique has been used extensively, despite the 
problem of establishing a proper Poincare surface of section 
in each case. However, the inspection becomes very difficult 
and also greatly deceiving in the case of dynamical systems 
with multidimensional phase space. 

(b) The maximal Lyapunov Characteristic Exponent (LCE) 
cr of an orbit informs us whether an orbit is ordered or chaotic. 
If cr > then the corresponding orbit is chaotic. Over thirty 
years ago, Benettin et al [2] studied theoretically the prob- 
lem of the computation of all LCEs and proposed an algo- 
rithm for their numerical computation. In particular, cr is 
computed as the limit for t —> oo of the quantity 



1 l|w(f)|| 
t l|w(0)||' 



(1) 



2 



Euaggelos E. Zotos 



where w(0) and w(f) are the deviation vectors for a given or- 
bit, at times / = and t > respectively. The time evolution 
of w is given by solving the so-called variational equations. 
Generally, for almost all choices of initial deviations w(0), 
the Umit for ? — > oo of equation (1) gives always the same 
value of (T 



cr — lim L, . 

/— >oo 



(2) 



In practice, of course, since the exponential growth of 
w(0 occurs for short time intervals, one stops the evolution 
of w(f) after some time Ti , records the computed Lj^ , nor- 
malize the vector w(?) and repeats the calculation for another 
time interval Ti, etc obtaining finally the cr as an average 
over many Ti, i= 1,2, N as 



1 " 

1=1 



(3) 



The basic drawback of the computation of cr is that, after 
every Ti, the calculation starts from the beginning and may 
yield an altogether different Lj. than the T(i-\) interval. Thus, 
since cr is influenced by the whole evolution of w(0), the 
time needed for Lt (or the L7 ) to converge is not known a 
priori and may become extremely long. This makes it often 
very difficult to determine whether cr finally tends to a pos- 
itive value (chaos) or converges to zero (order). The main 
advantage of the LCE is that it can be applied easily to dy- 
namical systems of any number of degrees of freedom. 

(c) The frequency analysis method proposed by Laskar 
[14-16, 19, 20], which is based on the calculation of the ba- 
sic frequencies of an orbit over a fixed interval of time. For 
orbits on the KAM tori, these frequencies are very accurate 
approximations of the actual frequencies of the dynamical 
system, but for chaotic orbits the computed values vary sig- 
nificantly in time and space. The frequency analysis method 
can be applied to systems with many degrees of freedom. 

(d) The study of the dynamical spectra of orbits [6, 21, 
28, 34, 35]. The distribution of the values of a given param- 
eter along the orbit has been proved a very reliable and pow- 
erful tool for the exploration of the properties of motion in 
Hamiltonian systems of two and three degrees of freedom. 
The reader can find more interesting information about the 
dynamical spectra in [35] and also in the references of this 
article. 

Over the last years, several methods have been intro- 
duced in order to characterize the nature of an orbit by study- 
ing the evolution of the deviation vectors, some of which 
are discussed in section 4. In the present paper we intro- 
duce and use a new, fast and easy to compute indicator: the 
Fast Norm Vector Indicator (FNVl). We focus our attention 
on the method of the FNVI, performing a systematic and 
thorough study of its behavior in the case of autonomous 
Hamiltonian systems with two (2D) and three (3D) degrees 
of freedom. The FNVI was found to fluctuate significantly 



for chaotic orbits, while it displays a nearly constant value 
for ordered orbits. 

This article is organized as follows: in section 2 we pro- 
vide the definition of the FNVI and we present results dis- 
tinguishing between ordered and chaotic motion in two and 
three degrees of freedom (2D and 3D) Hamiltonians, com- 
paring also the efficiency of the FNVI with the computation 
of the corresponding LCE. In the same section we establish 
a numerical criterion, that is the dFNVl, in order to quan- 
tify the results obtained by the FNVI method. In section 
3 we demonstrate the ability of the new method to reveal 
the detailed structure of the dynamics in the phase space in 
both 2D and 3D dynamical systems. In the following sec- 
tion, we apply the new method to follow the evolution of 
a sticky orbit in the two-dimensional dynamical system. In 
section 5 we conduct a detailed comparison between the new 
FNVl/dFNVl method and with some other relatively mod- 
em weU-known methods of chaos detection. FinaUy, in sec- 
tion 6 we summarize our results and we present a discussion 
and also the conclusions of the present research. 



2 The definition of tlie FNVI metliod 

The basic idea behind the FNVI method is the introduction 
of a simple and easily computed quantity that clearly iden- 
tifies the ordered or chaotic nature of an orbit. The FNVI is 
defined as 



FNVI(t) 



1 



|x(t)|| - ||x(0)|| 



l|x(0)|| 



t<t„ 



(4) 



where t is the time, || - || denotes the EucUdean norm of the 
vector x(t), while ||x(0)|| and ||x(f)|| are the norm vectors for 
a given orbit, at times f = and f > respectively. In prac- 
tice, we stop the evolution of x(f) after some time t\ - f, 
we record the computed FNVI(ti) and then we repeat the 
calculation for another time interval t2 = f, etc obtaining fi- 
nally the FNVI as a summation over many f,-, i - 1,2, A^. 
Here we must point out, that /* is the predefined time step of 
the numerical integration which remains constant during the 
entire predefined total integration time t„,ax- 

In order to apply the FNVI method, we shall investigate 
the nature of orbits in a dynamical system of perturbed har- 
monic osciUators given by the 3D potential 

Vix, y,z)=Y +y^+z'') + e {xY + y\'' + ^e-) , (5) 

where co is the common frequency of the oscillations along 
the X, y and z axis, while e is the strength of the perturbation. 
Potential (5) represents three coupled harmonic oscillators 
in the case of the 1:1:1 resonance. Potentials of this type are 
also known as perturbed elliptic osciUators [1, 3, 7]. The ba- 
sic reason for the choice of potential (5) is that perturbed 
elliptic oscillators appear very often in galactic dynamics 
and also in atomic-particle physics [7]. A second reason is 
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that is displays exact periodic orbits, interesting sticky or- 
bits together with large unified chaotic domains. Therefore, 
it gives us a great opportunity to test and prove the effective- 
ness and the reliability of the new FNVI method. 

The Hamiltonian function to the potential (5) reads 

H3 [x, y, z, Px,Py, Pz) =\{pI+pI+ pfj + Vix, y, z) 

= h, (6) 

where px, Py and p^ are the momenta per unit mass conju- 
gate to X, y and z respectively, while is the numerical of 
the Hamiltonian. 

The corresponding Hamiltonian for the 2D system can 
easily be obtained if we set z = = in equation (6). Then 

Hi [x,y,px,Py) =\{pI+ Py) + Vix,y) 

= ^{pI + Py + + (^y^) + ^^y^ 



2 

hi. 



(7) 



where hi is the numerical value of the 2D Hamiltonian. 

The outcomes of the present research are mainly based 
on the numerical integration of the equations of motion 

dyix,y,z) 



y = - 



z = -- 



dV(x,y,z) 



dz 



= -p^2e(x2+/)]z. 



(8) 



where the dot indicates derivative with respect to the time. 
We integrated numerically the equations of motion with a 
double precision Bulirsch-Stoer algorithm in FORTRAN 95 . 
The accuracy of our results was checked by the constancy of 
the energy integrals (6) and (7), which were conserved up to 
the fifteenth significant decimal point. 

We keep the involved parameters of the two systems 
fixed at the values 10 - \ and e = 1, while the value of 
the energy hi or h^ is treated as a parameter. 

A simple qualitative way of studying the dynamics of a 
Hamiltonian system is by plotting the successive intersec- 
tions of the orbits with a Poincare surface of section [17]. 
This method has been extensively applied to 2D Hamiltoni- 
ans, as in these systems the PSS is a two-dimensional plane. 
In 3D systems, however, the PSS is four dimensional and the 
behavior of the orbits cannot be easily visualized. One way 
to overcome this problem is to project the PSS to spaces 
with lower dimensions (see, [31, 32]). However, even these 
projections are often very complicated and difficult to be in- 
terpreted. 

In order to illustrate the behavior of the FNVI in 2D and 
3D dynamical systems, we first consider some representa- 
tive ordered and chaotic orbits. In Figure la we plot the in- 
tersection points of an ordered and a chaotic orbit of the dy- 
namical system (7), with a {x, px) PSS defined hyy = Q and 



Py > 0. The points of the ordered orbit lie on a torus and 
form a smooth closed curve on the PSS. On the other hand, 
the points of the chaotic orbit appear randomly scattered. 
The outermost black solid line shown in Fig. la is the Zero 
Velocity Curve (ZVC). The equation of the limiting curve 
ZVC (that it the curve containing all the invariant curves of 
the 2D system) is defined by the equation 

f2(x,Px)^^pl + V(x)^h2. (9) 

The time evolution of the FNVI for these two orbits and for 
a time period of lO' time units is plotted in Figure lb. In 
the case of the ordered orbit the FNVI remains almost con- 
stant, while in the case of the chaotic orbit it displays large 
fluctuations. The initial conditions for the regular orbit are: 
X() - 0.27, >'() = and pxo - 0, while the initial value of 77,0 
is found from the energy integral (7). The chaotic orbit has 
initial conditions: xq = 0.82, yo = and pxo = 0. For both 
orbits the value of the energy is hi = 1, while the the PSS 
shown in Fig. la was integrated for a time period of 5 x 10^ 
time units. In order to have a better and more closer view 
of the evolution of the FNVI for these two orbits, we present 
separately in Figure 2a-b the two diagrams. We observe, that 
in Fig. 2a which corresponds to the regular orbit, the FNVI 
after a small transient period of fluctuation for about 200 
time units, it stabilizes at a value and remains almost con- 
stant. On the contrary, in the case of the chaotic orbit shown 
in Fig. 2b the FNVI displays a completely different charac- 
ter with large and abrupt fluctuations throughout the entire 
period of the numerical integration. 

One may reasonably assume that the shape of the FNVI 
depends on the particular value of the total time of the nu- 
merical integration. Thus, if we integrate the FNVI of the 
chaotic orbit shown in Fig. 2b for a larger time period it 
might displays a similar pattern as the one of the ordered or- 
bit shown in Fig. 2a. In order to clarify this issue, we com- 
puted the FNVI for the same two 2D orbits discussed in Fig. 
2a-b but for a time period of 10^ time units. The results are 
presented in Figs. 3b and 3d respectively and compared with 
the corresponding LCE in order to test their validity. In Fig- 
ure 3a we see the evolution of the LCE of the ordered orbit. 
As expected the value of the LCE tends to zero as the time 
evolves. The evolution of the FNVI is presented in Figure 
3b. Again, the FNVI after a small transient period of fluc- 
tuation for about 200 time units, it stabilizes at a value and 
remains almost constant. This indicates the regular character 
of the orbit. On the other hand, the computation of the max- 
imal LCE, using Eqs. (1) and (2), despite its usefulness in 
many cases, does not have the same convergence properties 
over the same time interval. This becomes evident in Figure 
3 c where we plot the evolution of the LCE for the chaotic 
orbit. The computation of the LCE up to f = 10"* or even 
up to r - 10^ time units, still shows no clear evidence of 
convergence. Although Fig. 3c suggests that the orbit might 
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Fig. 1 (a-b): (a-left): The PSS of an ordered and a chaotic orbit of the Hamiltonian system (7). The ordered orbit corresponds to a closed (solid) 
elliptic curve, while the chaotic one is represented by the dots scattered over the PSS. The outermost black solid line is the ZVC. (b-right): The 
time evolution of the FN VI for the two orbits of panel (a). 
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Fig. 2 (a-b): The time evolution of the FNVl for a time interval of 10^ time units for a (a-left): regular 2D orbit and (b-right): chaotic 2D orbit. 
The initial conditions of the orbits and more details are given in the text. 



probably be chaotic since LCE remains different from zero 
for large time intervals, it does not allow us to conclude its 
chaotic nature with certainty, and so further computation of 
the LCE is needed. Thus, it becomes evident that an advan- 
tage of the FNVI, with respect to the computation of the 
LCE, is that the current value of the FNVI is sufficient to 
determine the chaotic nature of an orbit, in contrast to the 
maximal LCE, where the whole evolution of the deviation 
vector affects the computed value of LCE. Figure 3d depicts 
the evolution of the FNVI for the chaotic orbit. Therefore, it 
becomes clear that the shape of the FNVI does not depend 
on the time of the numerical integration. We observe in Fig. 
3d that the FNVI displays once more large and abrupt fluctu- 
ations throughout the time evolution. Note that in Figs. 3a-d 
the horizontal axis corresponding to the time is on log scale. 



in order to visuaUze more clearly the evolution of the indica- 
tors throughout the entire time interval. As the criterion by 
which we judge using the FNVI method, whether an orbit is 
regular or chaotic is quaUtative rather than quantitative, we 
could be sure for the chaotic nature of the orbit after only 
10^ time units of numerical integration. 

Let us now present some results about the 3D Hamil- 
tonian system (6). In Figure 4a we observe the evolution 
of the LCE of an ordered 3D orbit with initial conditions: 
xo = 0.09, Jo = 0,zo - 0.1,/?.ifl - p-s - 0, while the ini- 
tial value of PyQ is found from the energy integral (6). As 
expected, the value of the LCE tends to zero as the time 
evolves. The evolution of the FNVI is presented in Figure 
4b. The FNVI after a small transient period of fluctuation 
for about 250 time units, it stabilizes at a value and remains 
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Fig. 3 (a-d): The time evolution for a time interval of 10' time units of (a-upper left): the LCE of the regular 2D orbit, (b-upper right): the FNVI 
of the same regular orbit, (c-lower left): the LCE of the chaotic 2D orbit and (d-lower right): the FNVI for the same chaotic orbit. Note that the t 
axis is in log scale. 



almost constant, indicating the regular character of the orbit. 
In Figure 4c we see the evolution of the LCE of a chaotic 3D 
orbit. In this case the initial conditions of the chaotic orbit 
are: xq = 0.86, yo = 0,zo = 0.1, p^o = Pzo = 0. The LCE 
remains different from zero, for a time period of 10^ time 
units, which implies the chaotic nature of the orbit. The cor- 
responding evolution of the FNVI is presented in Figure 4d. 
It is clear that the pattem presented in Fig. 4d is completely 
different from that shown in Fig. 4b, where the 3D orbit is 
regular. For the 3D chaotic orbit, the FNVI displays large 
fluctuations and a highly irregular shape. Note that the t axis 
in Figs. 4a-d is in log scale. For both 3D orbits the value of 
the energy is /is = 1. 

Here, we have to point out that the criterion by which we 
determine so far the regular or chaotic nature of an orbit (2D 
or 3D), using the FNVI method is purely qualitative. In the 
case of a regular orbit the FNVI after a small transient pe- 
riod of fluctuation it stabilizes at a value and remains almost 
constant. On the other hand, when the orbit is chaotic the 
evolution of the corresponding FNVI follows a completely 
diff'erent pattem. It does not shows any evidence of conver- 
gence and displays large and abrupt fluctuations. This crite- 
rion has been obtained by testing a large number of regular 
and chaotic orbits in both dynamical systems (2D and 3D). 



Of course, since our criterion is qualitative we have to in- 
spect the shape of FNVI each time in order to characterize 
an orbit. Obviously, this is not very practical when someone 
wants to check a large volume of orbits, so as to form an idea 
about the global structure of the dynamical system. Thus, we 
need to establish a new numerical criterion, in order to quan- 
tify the results obtained by the FNVI method. This criterion 
can be derived by looking the shape of the FNVI shown in 
Figs. 3b, 3d, 4b and 4d. One may observe, that when the or- 
bit is regular the FNVI remains almost constant, while in the 
case of a chaotic orbit it displays high fluctuations. We are 
going to exploit this significant difference in order to obtain 
ovir quantitative criterion. Thus, we calculate the maximum 
and the minimum value of FNVI when t € [200, 1000]. We 
take this particular time interval because in the case where 
the orbit is regular, for t < 200 time units there is a transient 
period of fluctuation, which may create a malfunction to our 
criterion. Using the above procedure we can define the 



dFNVI = FNVIn,ax - FNVI„ 



(10) 



where FNVImax and FNVImin are the maximum and the min- 
imum value of FNVI respectively when t e [200, 1000]. The 
value of dFNVI can provide us the quantitative criterion that 
we seek. We note, that it is not easy to define a threshold 
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Fig. 4 (a-d): The time evolution for a time interval of 10' time units of (a-upper left): the LCE of a regular 3D orbit, (b-upper right): the FNVI of 
the same regular orbit, (c-lower left): the LCE of a chaotic 3D orbit and (d-lower right): the FNVI for the same chaotic orbit. Note that the t axis 
is in log scale. The initial conditions of the orbits and more details are given in the text. 



value, so that the dFNVI being larger than this value reli- 
ably signifies chaoticity. Nevertheless, extensive numerical 
experiments in both dynamical systems indicate that in gen- 
eral, a good guess for this value could be 0.05. Thus, when 
dFNVI > 0.05 the orbit is chaotic, while when dFNVI < 0.05 
the orbit ir ordered. This threshold value appUes in both 2D 
and 3D orbits. For the 2D regular orbit presented in Fig. 2a 
dFNVI = 0.009, while for the 2D chaotic orbit of Fig. 2b 
dFNVI - 0.45. Similarly, for the 3D regular orbit presented 
in Fig. 4b dFNVI = 0.006, while for the 3D chaotic orbit 
of Fig. 4d dFNVI = 0.52. Therefore, we may conclude that 
the quantification of the FNVI method has been archived 
successfully. In the next section, we shall provide more ex- 
tensive and detailed results using the dFNVI. 



3 Distinguishing between regions of order and chaos 

The dFNVI offers indeed a very easy and efficient method 
for distinguishing the chaotic versus ordered nature of or- 
bits in a variety of problems. In the present section, we use 
it for identifying regions of phase space where large scale or- 
dered and chaotic motion are both present. In particular, we 
shall apply the dFNVI method in order to reveal the struc- 



ture in both Hamiltonian systems of two (2D) and three (3D) 
degrees of freedom. 



3.1 Order and chaos in the 2D dynamical system 

In Figure 5a-d we present detailed plots of the (x, px), y -0, 
Py > PSS of the 2D dynamical system (7) for different 
values of the energy h2. Fig. 5a shows the PSS plot, when 
/z2 = 1. One can see, that regions of ordered motion around 
stable periodic orbits are seen to coexist with chaotic re- 
gions filled by scattered points. In particular, a large part of 
the phase plane is covered by chaotic orbits, while the reg- 
ular regions are occupied by invariant curves mainly around 
the points ixo,Pxo) = (0,0) and (xo,/?io) - (0,±-\/h2). The 
above points give the position of two exact periodic orbits. 
The first is the y axis, while the second describes thex= ±y 
straight-line orbits on the (x, p^) phase plane. Furthermore, 
one observes a large number of smaller islands of invariant 
curves produced by secondary resonances and some sticky 
regions as well. Fig. 5b is similar to Fig. 5a but when h2 = 2. 
Here, we see that the majority of the phase plane is covered 
by a unified chaotic sea. The regular regions are confined 
only around the straight line periodic orbits, while the y axis 
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has now become unstable. Some smaller islands of invari- 
ant curves are also present. The most interesting feature ob- 
served in this case, is the sticky regions around the chain of 
islands of invariant curves on the x as well as on the px axis. 
If we increase the value of energy to h2 = 4, we obtain the 
phase plane shown in Fig. 5c. In this case, the >' axis has 
returned to stabihty, while the straight line periodic orbits 
are unstable. Almost all the phase plane is chaotic, except of 
a small region around the center and some small islands of 
invariant curves which are products of secondary resonant 



orbits. A sticky region around the chain of the small islands 
of invariant curves near the center is also observed. Fig. 5d 
is similar to Fig. 5a, but when h2 = 7.2. Here, almost the 
entire phase plane is covered by chaotic orbits. There is a 
very small regular region confined near the center. A care- 
ful observation also reveals some very tiny regular islands 
of invariant curves embedded in the vast chaotic sea. 

In order to demonstrate and prove the effectiveness of 
the dFNVI method, we first consider orbits whose initial 
conditions he on the Unes x = and = 0. In particular, 
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Fig. 6 (a-h): The values of the dFNVI for orbits of the 2D system (7) with initial conditions on the (left pattern): p,- = line as a function of the 
xq coordinate of the initial conditions (xq, 0), (right pattern): x = line as a function of the p,o coordinate of the initial conditions (0, Pxo), for the 
PSS plots of Fig. 5a-d. Green and red dots indicate initial conditions corresponding to regular and chaotic orbits respectively. See text for more 
details. 
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we take 5 x lO-' equally spaced initial conditions on these 
lines and we compute the value of the dFNVI for each one. 
The results are presented in Figure 6a-h where we plot the 
dFNVI as a function of the xq and /?,o coordinates of the ini- 
tial conditions of these orbits for t = 10^ time units. In all 
panels the data points are line connected, so that the changes 
of the dFNVI values are clearly visible. With green dots we 
represent the initial conditions (xq, 0) or (0, p^o) correspond- 
ing to regular orbits, while the initial conditions producing 
chaotic orbits are marked with red dots. Note that there are 
intervals where the dFNVI has low values (e.g. lower than 
0.05), which correspond to ordered motion inside the islands 
of stability crossed by the x = and Px = lines shown 
in Fig. 5a-d. There also exist regions where the dFNVI has 
large values (e.g. larger than 0.1) denoting that in these re- 
gions the motion is chaotic. These intervals correspond to 
the regions of scattered points crossed by the x = and 
Px = lines in Fig. 5a-d. Although most of the initial condi- 
tions give large (> 0.1) or very small (< 10"^) values for the 
dFNVI, there also exist initial conditions that have interme- 
diate values of the dFNVI (0.05 < dFNVI < 0.1). These ini- 
tial conditions correspond to sticky chaotic orbits, remain- 
ing for long time intervals at the borders of islands, whose 
chaotic nature will be revealed later on. In Figs. 6a and 6b 
we observe the values of the dFNVI computed for t = 10^ 
time units for orbits with initial conditions on the Px = and 
X = lines on the PSS shown in Fig. 5a, as a function of the 
xo coordinate and the pxo momenta respectively. Similarly, 
Figs. 6c and 6d correspond to the PSS of Fig. 5b, Figs. 6e 
and 6f correspond to the PSS of Fig. 5c and Figs. 6g and 6h 
correspond to the PSS shown in Fig. 5d. 

In Fig. 6a, around xq ~ 1.05 there exists a group of 
points inside a large chaotic region having dFNVI < 0.05. 
These points correspond to orbits with initial conditions in- 
side a small stabiUty island, which is not very visible in the 
detailed PSS plot of Fig. 5a. Also the points with pxo = ± 
0.91 in Fig. 6b, have been characterized by the dFNVI as 
initial conditions corresponding to chaotic orbits, while all 
their neighboring points have dFNVI < 0.05. These points 
actually correspond to a weak chaotic separatrix inside the 
vast domain of stability, which can be revealed only after a 
very high magnification of this region of the PSS. So, we see 
that the systematic application of the dFNVI method can re- 
veal very fine details of the structure in a dynamical system. 

It will be of particular interest, to calculate the mean 
value of dFNVI for the regular and chaotic orbits discussed 
in Fig. 6a-h. In other words, we are going to calculate the 
mean value of dFNVI for the regular orbits, that is (dFNVI)R 
and the mean value of dFNVI for the chaotic orbits (dFN VI)c, 
with initial conditions along the x and the Px axis for differ- 
ent values of the energy h2- In Table 1 we present our re- 
sults. It is evident, that as the value of the energy increases, 
we have an increase at the mean value of the dFNVI for 



both regular and chaotic orbits. In the case of chaotic orbits, 
the increase of (dFNVI)c indicates that the chaoticity of the 
dynamical system increases as the value of the energy is am- 
plified. We shall come back to this point later in this section. 



Table 1 The mean value of dFNVI for the regular and chaotic orbits 
discussed in Fig. 6a-h. 





X - axis 


Px- 


axis 


hi 


<dFNVI)R 


(dFNVI)c 


(dFNVI>R 


(dFNVI)c 


1.0 


0.00725 


0.32251 


0.00551 


0.28687 


2.0 


0.01136 


1.14872 


0.01922 


1.15079 


4.0 


0.01732 


1.27763 


0.02317 


1.21792 


7.2 


0.02181 


1.43602 


0.02761 


1.42706 



By carrying out the previously presented analysis for 
initial conditions (xo,Pxo) not only along a line but on the 
whole phase plane of a PSS and giving to each point a color 
according to the value of the dFNVI, we can have a clear 
picture of the regions where chaotic or ordered motion oc- 
curs. The outcomes of this procedure for the 2D dynamical 
system (7), using a dense grid of initial conditions (xcp^o) 
on the PSS, are presented in Figure 7a-d. The value of the 
energy h2 in Figs. 7a, 7b, 7c and 7d is the same as in Figs. 5a, 
5b, 5c and 5d respectively. Thus, in Figs. 7a-d we clearly dis- 
tinguish between green regions, where the motion is ordered 
and red regions, where it is chaotic. The outermost black 
soUd Une shown in Figs. 7a-d is the hmiting curve (ZVC) 
defined by Eq. (9). It is worth mentioning, that in Fig. 7a we 
can observe small islands of stabihty inside the large chaotic 
sea, which are not very visible in the detailed PSS of Fig. 5a, 
such as that for xq ~ + 1.05, pxo ~ 0. Although all the or- 
bits with initial conditions (xo,/?,to) in Figs. 7a-d were com- 
puted for only / = 10^ time units (such as Figs. 6a-h), this 
time was suflicient enough for the clear revelation of small 
ordered regions inside the chaotic domains. We must point 
out, that the outcomes derived using the dFNVI method re- 
garding the structure of the phase plane coincide with those 
obtained using the PSS technique. In other words, we see 
that the structure of the phase planes shown in Figs. 7a-d are 
practically identical with those presented in Figs. 5a-d. For a 
grid of about 10^ equally spaced initial conditions (xo,pxo), 
we need about 7 h of CPU time on a Pentium Dual Core 
Processor at 2.2GHz PC, in order to construct each grid-plot 
shown in Figs. 7a-d. 

Grid-plots such as those of Fig. 7a-d, apart from pre- 
senting the regions of order and chaos, can also be used to 
provide very accurate estimations regarding the fraction of 
the phase-space volume occupied by chaotic or ordered or- 
bits and also give us good initial guesses for the location 
of stable periodic orbits, in regions where the motion is or- 
dered. One can conclude form Figs. 5a-d or from Figs. 7a-d, 
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Fig. 7 (a-d): Regions of different values of the dFNVI on the PSS plots of the 2D dynamical system (7) when (a-upper left): hi = 1, (b-upper 
right): h2 = 2, (c-lower left): h2 = 4 and (d-lower right): hi = 7.2. In all panels, the initial conditions (xq, Pm) are colored green if their dFNVI < 
0.05 and red if dFNVI > 0.05. 



that the fraction of the phase plane covered by chaotic orbits 
increases as the value of the energy h2 increases. Figure 8 
shows a plot of the percentage of the phase plane A% cov- 
ered by chaotic orbits versus h2. We observe that the per- 
centage A% increases rapidly, as the value of h2 increases. 
Dots indicate the values of the percentage A% obtained nu- 
merically from the grid-plots, while the solid line is a fourth 
degree polynomial fitting curve. We must point out that the 
percentage A% is calculated as follows: when constructing 
the grid-plots, we count the initial conditions (-«o,/'.rf)) cor- 
responding to regular orbits and also the initial conditions 



(xo,Px-o) corresponding to chaotic orbits. Then, we divide 
the number of chaotic orbits by the total number of orbits 
and thus, we obtain the percentage A% for each value of 
the energy h2. Note, that in every case the total number of 
the calculated orbits is about 10^, while the exact number of 
chaotic orbits varies depending on the value of the energy 

h2. 

In order to have an estimation of the degree of chaos of 
the dynamical system from another point of view, we have 
plotted the average value of the maximal LCE versus h2. 
The results are shown in Figure 9a. Note, that the ( LCE ) 
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Fig. 9 (a-d): (a-left): A plot of the average value of the maximal LCE versus A2 and (b-right): A plot of the < dPNVI > as a function of the value of 
the energy /12. 
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Fig. 8 A plot of the percentage of the phase plane A% covered by 
chaotic orbits versus the value of the energy h2. 

increases linearly with ^2- Here, we must point out that it 
is well known that the value of the LCE is different in each 
chaotic component [22]. As we have in all cases regular re- 
gions and only one unified chaotic sea in each (x, px) phase 
plane, we calculate the average value of LCE by taking 10^ 
orbits with different and random initial conditions (xo,/?a()) 
in the chaotic domain for every value of the energy h2. Note 
that all calculated LCEs are different on the fifth decimal 
point in the same chaotic region. We follow the same pro- 
cedure using now the dFNVI. For the same 10^ initial con- 
ditions corresponding to chaotic orbits for each value of the 
energy /z2, we calculated the ( dFNVI ). The results are pre- 
sented in Figure 9b. Again the ( dFNVI ) increases linearly 
with h2. The increase of { dFNVI ) indicates that the chaotic- 
ity of the dynamical system increases as the value of the en- 
ergy is amplified. As the behavior of ( dFNVI ) coincides 
with the behavior of ( LCE ), we prove once more the reli- 
abiUty of dFNVI method. However, the main advantage of 
the dFNVI method is that it requires only 10^ integration 
time units, while the LCE needs at least about 5 x 10"* in- 
tegration time imits so as to provide reUable and definitive 



evidence. Thus, the dFNVI is a much more faster indicator 
than the LCE. 

3.2 Order and chaos in the 3D dynamical system 

In the case of 3D Hamiltonian systems the PSS is now four 
dimensional and thus, not so useful as in the 2D systems. On 
the other hand, the dFNVI method has the ability once more 
to identify successfully regions of order and chaos in the 
phase space. To prove this, let us start with initial conditions 
on a 4D grid of the PSS. In this way, we find again regions of 
order and chaos, which may be visuaUzed, if we restrict our 
study to a subspace of the whole 6D phase space. We con- 
sider orbits with initial conditions (xq, Zo.Pja). Jo = /'zO = 0> 
while the initial value of pyo is always obtained from the en- 
ergy integral (6). In particular, we define a value of zq which 
is kept constant and then we calculate the dFNVI of 3D or- 
bits with initial conditions ixo,pxa), yo = PzO - 0. Thus, we 
are able to construct again a 2D plot depicting the (x, px) 
plane but with an additional value of zo- All the initial con- 
ditions of the 3D orbits lie inside the limiting curve defined 
by 

1 , 

/3(x,p.,;z,)) = i^Px + V(x;zo) = hj,. (11) 

Our results are presented in Figure 10 a-d, where we give 
four grid-plots for the same value of the energy = \, but 
for different initial value of zo- Again, we can distinguish 
between regions of ordered (colored in green, where dFNVI 
< 0.05) and chaotic motion (colored in red, where dFNVI > 
0.05). In Fig. 10a where zo - 0.01 we see that the structure 
of the (x, px) plane is quite similar with the corresponding 
2D grid-plot shown in Fig. 7a. We also see that well defined 
islands of stabiUty inside the unified chaotic sea still exist. 
The main difference from the 2D grid-plot of Fig. 7a is that 
now in the 3D system, a large portion of orbits around the 
stable periodic points on the px axis have altered their nature 
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Fig. 10 (a-d): Regions of different values of the dFNVI of the 3D dynamical system (6) for = I, when (a-upper left): zo = 0.01, (b-upper right): 
Zo = 0.1, (c-lower left): zo = 0.3 and (d-lower right): zo = 0.5. In all panels, the initial conditions {xo,zo, Pxo) are colored green if their dFNVI < 
0.05 and red if dFNVI > 0.05. 



form ordered to chaotic. Fig. 10b is similar to Fig. 10a but 
when Zo = 0.1. In this case, the structure of the (x, px) plane 
is very similar to that shown in Fig. 10a. One may observe, 
that the percentage corresponds to chaotic orbits has been 
increased and moreover the islands of stability have begun 
to destabilize and lose their well defined structure. In Fig. 
10c we increase the initial value of zo to 0.3. In this case, 
it is evident that the initial conditions correspond to ordered 
orbits are delocalized, since now we can see no signs of well 
defined islands of stability. Finally, in Fig. lOd we present a 
grid-plot when zq - 0.5. Here, the vast majority of the ini- 



tial conditions correspond to chaotic orbits, while the initial 
conditions correspond to ordered orbits are completely de- 
localized and randomly scattered mainly in the central part 
of the {x, Pjc) plane. All the 3D orbits with initial conditions 
{xo,zo,Pjco) in Figs. lOa-d were computed for only 10"^ time 
units. With a closer look to Figs. lOa-d one may conclude 
the following two main points: (i) the increase of the ini- 
tial value of Zo has as a result the increase of the percent- 
age on the (x, p^) plane corresponding to chaotic orbits and 
(ii) as we amplify the value of zo we observe that the entire 
area of the (x,px) plane defined by Eq. (11) is reduced and 
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Fig. 11 A plot depicting the relationship between the percentage A% 
of the initial conditions (xo,pm) corresponding to chaotic 3D orbits 
and the initial value of zo, when hs = 1. 

becomes smaller. This can easily be explained. The coordi- 
nates (x, y, z) and the momenta (p^, Py, Pz) of the 3D orbits 
are connected through the energy integral (6) and since we 
increase the value of zo while keeping constant the value of 
the energy this has as the result the permissible values 
of the initial conditions ixo,Pxo) to be reduced. For a grid 
of about 10'' equally spaced initial conditions (xo,/?xo), we 
need about 7.5 h of CPU time on a Pentium Dual Core Pro- 
cessor at 2.2GHz PC, in order to construct each grid-plot 
shown in Figs. lOa-d. 

The grid-plots shown in Figs. lOa-d could be considered 
as (x, px) "phase planes" but not with the strict definition. In 
fact, these plots show the structure of the (x, px) subspace 
for a given value of the zo- Therefore, we study the subspace 
(x,px), z = Zo, y = Pz = with > of the 6D phase 
space of the 3D dynamical system. Note that grid-plots such 
as those of Figs. lOa-d, apart from presenting the initial con- 
ditions correspond to ordered and chaotic orbits, can also be 
used to estimate the fraction of the phase-space volume oc- 
cupied by chaotic or ordered 3D orbits. In Figure 11, we 
present a plot depicting the relationship between the per- 
centage A% of initial conditions {xo,Pxo) corresponding to 
chaotic 3D orbits and the initial value of zo- We observe that 
the chaotic percentage zo increases linearly with zo- Our nu- 
merical calculation suggest, that when zo ^ 0.67 for = 1 
eventually all the initial conditions ixo,zo,Pxo) correspond 
to chaotic orbits, while regular motion, if any, is negligible. 
Note, that the plot shown in Fig. II corresponds to = I. 
For different values of the energy the relationship between 
A% and zq is quite similar, but as we amplify the value of 
the energy h^, the increase of the chaotic percentage is more 
rapid and does not follow a linear pattern any longer. 

Finally, in order to have an estimation of the degree of 
chaos of the 3D dynamical system from another point of 
view, we have plotted the average value of the maximal LCE 
versus the value of the energy h^. The results are shown in 



Figure 12a. Note that the ( LCE ) increases linearly with h^, 
when the initial value of zo is O.OI. As we have regular re- 
gions and only one unified chaotic sea in the (x, Px) planes 
(see Fig. 10a), we calculate the average value of the LCE 
by taking 10^ orbits with difl'erent and random initial con- 
ditions (xq,Pxq,zq - 0.01), >'o = Pzo - and pyo > in 
the chaotic domain for every value of the energy hj,. Note, 
that all calculated LCEs are different on the fifth decimal 
point in the same chaotic region. For different initial value 
of the Zo, our numerical experiments indicate that the rela- 
tionship between ( LCE ) and remains linear and similar 
to that shown in Fig. I2a. Here we should point out, that the 
{ LCE ) of the 3D system is smaller than the corresponding 
( LCE ) of the 2D system shown in Fig. 9a. We follow the 
same procedure using now the dFNVI method. For the same 
10^ initial conditions ixo,Pxo,zo = 0.01), yo = Pzo = and 
PyO > corresponding to chaotic 3D orbits for each value 
of the energy /i3, we calculate the ( dFNVI ). The results are 
presented in Figure 12b. Again, the < dFNVI ) increases lin- 
early with hi. The increase of the ( dFNVI ) indicates that 
the chaoticity of the dynamical system increases as the value 
of the energy is amplified. For different initial value of the 
Zo, our numerical calculations indicate that the relationship 
between ( dFNVI ) and hj, also remains linear and similar 
to that shown in Fig. 12b. As the behavior of ( dFNVI ) 
coincides with the behavior of { LCE ), we prove the reU- 
ability of dFNVI method also in the 3D dynamical system. 
Of particular interest is the fact, that the < dFNVI ) of the 3D 
system is also smaller than the corresponding { dFNVI ) of 
the 2D system shown in Fig. 9b. This strongly suggests that 
the chaoticity in the dynamical system of three degrees of 
freedom(3D) is smaller that in the corresponding dynamical 
system of two degrees of freedom (2D), when h2 = h^. 



4 Application in the case of a sticky orbit 

In this section, we shall try to test the FNVI method's ability 
posing a somehow more difficult task, that is to distinguish 
between ordered and chaotic motion, in the case of chaotic 
orbits that seem to be regular for a number of periods. Such 
orbits, called "sticky orbits" and they usually grow near the 
outer regions of an island of stabiUty [18]. Cantori that ex- 
ist outside that region may restrict the motion in a very thin 
layer for a long time interval. After that time interval, the or- 
bit can escape to the surrounding chaotic region, through the 
gaps of the cantorus. Due to this peculiar behavior, a sticky 
orbit often resembles an ordered one seen in a Poincare sur- 
face of section, until the orbit escapes to the chaotic domain. 

We follow the evolution of a two-dimensional sticky or- 
bit which produces two sets of five small islands of invari- 
ant curves near the x axis, using the PSS technique. This 
sticky orbit has initial conditions: xq = -1.705, yo = and 
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Fig. 12 (a-d): (a-left): A plot of the average value of the maximal LCE 
of the energy hj, when zo = 0.01. 



= 0, while the value of Pyo is obtained from the en- 
ergy integral (2). Our results are presented in Fig. 13 (a-d). 
Fig. 13a shows the regions formed in the (x, p^) phase plane 
by this orbit, for a period of 1200 time units of numerical 
integration. We see two sets of five separate islands of in- 
variant curves. The outermost black solid line shown in Fig. 
13a is the Zero Velocity Curve (ZVC) defined by Eq. (9). 
When T - 1200 time units (Fig. 13a), one can see that the 
orbit is chaotic merely zooming in the plot and seeing that 
the points are not smoothly distributed over a well defined 
one-dimensional curve. In Fig. 13b the same orbit was cal- 
culated for 1700 time units. Now, we observe that in each 
case the five small islands of invariant curves were joined 
together. By observing Figs. 13a and 13b, it is obvious that 
from the pictures presented in the PSS no safe conclusion 
can be reached concerning the regular or chaotic behavior 
of this orbit. Fig. 13c depicts the same orbit calculated for 
a time period of 12500 time units. It is evident, that now 
the test particle has left the sticky region and begins to en- 
ter the vast surrounding chaotic region. In Fig. 13d we have 
calculated the orbit for 5x10"* time units. The chaotic na- 
ture of the orbit is now completely revealed, since the scat- 
tered points in the {x, px) phase plane fill the entire available 
chaotic domain. 

A better and more enlightening view for the evolution of 
the sticky orbit can be seen using the S (g) dynamical spec- 
trum [33,35]. Fig. 14a shows the S (g) spectrum of the sticky 
orbit for 1 150 time units of numerical integration. Here we 
observe ten separate t/-type structure with a large number 
of additional small and large peaks indicating the sticky na- 
ture of the orbit. In Fig. 14b we present the S(g) spectrum 
computed for a time interval of 1670 time units. In this case, 
the ten structures have joined together producing two sep- 
arate new structures. Fig. 14c depicts the S (g) spectrum of 
the same orbit calculated for 12200 time units of numerical 
integration. One may observe, that the two separate struc- 
tures shown in Fig. 14b have now joined together producing 
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a single spectrum. The shape of this spectrum has all the 
characteristics of a chaotic spectrum and thus strongly in- 
dicates that the test particle has undoubtedly left from the 
sticky region and continued its journey to the surrounding 
chaotic sea. Fig. 14d shows the S(g) spectrum computed for 
5x10'' time units. We see that the spectrum has the usual 
shape of a spectrum corresponding to a chaotic orbit. Here 
we must point out, that the time intervals regarding the grad- 
ual time evolution of the sticky orbit obtained by the S (g) 
spectrum are very close to those derived by the formation of 
the regions in the {x, p^) phase plane, shown in Fig. 13 (a-d). 



Now let's use some other dynamical indicators in order 
to follow the evolution of the sticky orbit and also compare 
and verify our results. Fig. 15 (a-b) shows the time evolu- 
tion of the LCE and the FNVI for the same sticky orbit dis- 
cussed and studied earlier. Fig. 15a shows the time evolution 
of the LCE. It is evident that after 12600 time units the orbit 
gradually becomes chaotic. In Fig. 15b we present the time 
evolution of the FNVI for this orbit. We observe, that the 
pattern of the time evolution of the FNVI is almost identical 
with the pattern of the LCE shown in Fig. 15a. Furthermore, 
the FNVI clearly indicates that after 12100 time units the 
orbit gradually alters its nature to chaotic. Note that in Fig. 
15a-b the t axis is in log scale in order to have a more clear 
view of the time evolution. Therefore, we may conclude that 
the results obtained by the PSS technique, the S (g) dynam- 
ical spectrum, the LCE and also the FNVI method shown 
in Figs. 13a-d, 14a-d, 15a and 15b respectively, coincide 
that the sticky orbit becomes chaotic after a time interval 
of about 12500 time units of numerical integration. Thus, 
we may conclude that the new FNVI method has the abil- 
ity to follow the time evolution of a two-dimensional sticky 
orbit and provide rehable results regarding the time needed 
for the transition from sticky motion to chaos. 



Detecting order and chaos in Hamiltonian systems by tlie Fast Norm Vector Indicator (FNVl) method 



15 




5 Comparison with other dynamical indicators 

The results presented in the previous sections reveal, that 
the dFNVI is a simple, efficient and easy to compute tool 
for the distinguishing between ordered and chaotic motion 
in Hamiltonian systems. Implementing the dFNVI method 
is a very easy computational task, as we only have to follow 
the evolution of the norm of an orbits 's vector throughout 
the integration time interval and compute the maximum and 
the minimum value of FNVI when t e [200, 1000]. In the 
case of regular motion the FNVI after a small transient pe- 
riod of fluctuation it settles down to a value and remains al- 



most constant, while the dFNVI is always smaller than 0.05. 
In particular, our extensive numerical calculations suggest, 
that when an orbits is regular the dFNVI is much smaller 
that 0.01. However, we decided to increase the threshold 
value to 0.05 in order to obtain more secure and reUable re- 
sults. On the other hand, in the case of chaotic motion the 
FNVI continues to fluctuate randomly without giving any 
sign of convergence, while the numerical criterion, that is 
the dFNVI is always larger than 0.05. It is exactly this dif- 
ferent behavior of the FNVI regarding the convergence, that 
makes it an ideal tool of chaos detection. However, the re- 
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Fig. 14 (a-d): Time evolution of the 5(g) dynamical spectrum for the sticky orbit. The integration time is (a-upper left): T = 1150 time units, 
(b-upper right): T = 1670 time units, (c-lower left): T = 12200 time units and (d-lower right): 5 x 10* time tinits. 



0.20 



0.15 



Q 0.10 



0.05 



0.00 




1 



4.5 



4.0 



sr 3.5 



3.0 



2.5 




1 



2 3 4 5 

a log(t) b 

Fig. 15 (a-b): Time evolution of the (a-left): the LCE and (b-right): the FNVI for the same sticky orbit. 
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suits obtained by the FNVI are only qualitative. Thus, we 
established a numerical criterion the dFNVI, which can pro- 
vide quantitative results about the nature of an orbit. As we 
have seen in the previous sections, the results obtained by 
the dFNVI are completely rehable, conclusive and beyond 
any doubt. The dFNVI helps us decide the chaotic or regu- 
lar nature of orbits faster and with less computational effort 
than the estimation of the maximal LCE. This happens be- 
cause the time needed for the LCE in order to give a clear 
and undoubted indication of convergence to non-zero values 



is usually much greater than the time needed for the dFNVI 
to provide reliable evidence regarding the character of an 
orbit. 

Many other dynamical indicators have been introduced 
in the recent years, some of which are compared in this sec- 
tion with the dFNVI method. In order to verify and prove, 
once more, the effectiveness and the reliabiUty of our new 
method, we shall compare our results with four other well 
known dynamical indicators: (a) The Fast Lyapunov Indi- 
cator (FLI), (b) The Relative Lyapunov Indicator (RLI), (c) 
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The Smaller Alignment Index (SALl) and (d) The General- 
ized Ahgnment Index (GALI). At this point we believe, that 
it should be wise to recall the definitions of these indicators, 
in order to help the reader to compare the results of each 
method. 

The FLI introduced in [8, 9] has been proved a very fast, 
reliable and efl'ective tool, which can be defined as 



FLI(t) = log||w(t)||, t<t„ 



(12) 



where vf(t) is a deviation vector. The computation of the FLI 
on a relatively short time t^ax is enough to discriminate be- 
tween chaotic and regular orbits. The FLI of a regular orbit 
increases linearly, while for a chaotic orbit, the FLI increases 
super-exponentially. Moreover, the FLI may be used to dis- 
criminate among regular motion between non resonant and 
resonant orbits. In order to avoid the linear trend of the FLI 
along regular orbits, we can use the de-trended FLI (DFLI), 
which is simply the FLI divided by t. Thus, 

l|w(t)|p 



DFLI(t) = log 



t<t„ 



(13) 



In this case, the classification of an orbit is based on the 
clearly distinct behavior of the DFLI pattern. In particular, 
in the case of a regular orbit the DFLI is almost constant 
with values lower than 10, while in the case of a chaotic 
orbit, we observe a very rapid increase of the DFLI. We be- 
lieve, that using definition (13) instead of (12) we can obtain 
better and more conclusive results. The DFLI is a relatively 
fast indicator, as it needs only about 5x10^ times units in 
order to provide reliable results regarding the nature of an 
orbit. Note, that for the calculation of the DFLI, we need 
the assistance of the variational equations and of course the 
computation of the deviation vector w(f). 

Another interesting dynamical indicator is the RLI pro- 
posed by Sandor et al [23]. The RLI is practically the ab- 
solute value of the difference of Li of two initially nearby 
orbits and can be defined as 



RLl(t) = |L (x H- 6x, t) - L (x, t)| , 



(14) 



where 



. . 1, Ilw2(r)ll 



Using the variational equations and the evolution of the 
deviation vectors we can compute the SALI [25, 26]. The 
basic idea behind the SALI method is the introduction of 
a simple quantity that clearly indicates if a deviation vec- 
tor is aligned with the direction of the eigenvector which 
corresponds to the maximal LCE. In general, any two ran- 
domly chosen initial deviation vectors wi(0) and W2(0) will 
become aligned with the most unstable direction and the an- 
gle between them will rapidly tend to zero [2]. Thus, we 
check if the two vectors have the same direction in phase 
space, which is equivalent to the computation of the above- 
mentioned angle. More specifically, we follow simultane- 
ously the time evolution of an orbit with initial condition 
x(0) and two deviation vectors with initial conditions wi(0) 
and W2(0). As we are only interested in the directions of 
these two vectors we normaUze them, at every time step, 
keeping their norm equal to 1 . This controls the exponen- 
tial increase of the norm of the vectors and avoids overflow 
problems. Since, in the case of chaotic orbits the normal- 
ized vectors point to the same direction and become equal 
or opposite in sign, the minimum of the norms of their sum 
iantiparallel alignment index) or difference {parallel align- 
ment index) tends to zero. So, the SALI is defined as 



SALI(t) = min(d_(t),d+(t)), 



where d-{t) - 



while \6iL\ 1. Very smaU values of the RLI (rLI < 10"") 

denote ordered motion, while larger values (rLI > lO '*) de- 
note chaotic behavior (for more information on the RLI method 
see [23]). The RLI method needs at least about 5x10'' times 
units in order to provide rehable results regarding the nature 
of an orbit. We must point out, that the computation of the 
RLI requires the time evolution of two orbits and two de- 
viation vectors {y/\{t) and y/2{f) one for each orbit), while 
the computation of the dFNVl is much faster as we compute 
only the evolution of the main orbit without the use of any 
deviation vector. 



W2(0 



and d+{t) - 



(16) 

W2(0 II 

|l|Wl(OII l|W2(t)ll H " + ||l|Wl(OII ' IIW2WIIII 

are the parallel and the antiparallel alignment indices respec- 
tively. From the above definition it is evident that SALl(t) 
e [0, V2] and when SALI = the two normalized vectors 
have the same direction, being equal or opposite. Imple- 
menting the SALI method is an easy computational task, as 
we only have to follow the evolution of an orbit and of two 
deviation vectors, computing in every time step the mini- 
mum norm of the difference and the addition of these vec- 
tors. In the case of chaotic motion the SALI eventually tends 
exponentially to zero, reaching rapidly very small values or 
even the limit of the accuracy of the computer (10"^^). On 
the other hand, in the case of ordered motion the SALI fluc- 
tuates around non-zero values. The SALI has a clear phys- 
ical meaning as zero, or a very small value of the index, 
signifies the alignment of the two deviation vectors. An ad- 
vantage of the method is that the index ranges in a defined 
interval (SALI e [0, "V^]) and so very small values of the 
SALI (e.g., smaller than 10"^) establish the chaotic nature 
of an orbit beyond any doubt. The SALI is a very fast indi- 
cator, as it needs only about 10^ times units in order to reveal 
the regular or chaotic character of an orbit. 

Recently, Skokos et al [27] generalized and improved 
considerably the SALI method mentioned above, by intro- 
ducing the GALI method. This indicator retains the advan- 
tages of the SALI i.e. its simplicity and efficiency in dis- 
tinguishing between regular and chaotic motion but, in ad- 
dition, is faster than the SALI, displays power law decays 
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that depend on torus dimensionality and can also be applied 
successfully to cases where the SALI is inconclusive, like in 
the case of chaotic orbits whose two largest Lyapunov ex- 
ponents are equal or almost equal. For the computation of 
the GALI we use information from the evolution of more 
than two deviation vectors with respect to the reference or- 
bit, while salt's computation requires the evolution of only 
two such vectors. In particular, GALIk is proportional to 
"volume" elements formed by k initially linearly indepen- 
dent unit deviation vectors whose magnitude is normalized 
to unity at every time step. If the orbit is chaotic, GALIt 
goes to zero exponentially fast by the law 

GALIk(t) oc e-K'^i-^2)+(<ri-cr3)+...+(<r,-(rt)]t_ ^^jy 

If, on the other hand, the orbit lies in an A^-dimensional 
torus, GALIk displays the following behaviors: Either 

GALIk(t) « constant for 2 < k < N, (18) 

or, if A'^ < A: < 2N, it decays with different power laws, de- 
pending on the number m of deviation vectors which initially 
lie in the tangent space of the torus, i.e., 
1 



GALIk(t) oc 



1 



if N<k<2N and 0<m<k-N 



if N <k<2N and m>k-N. 



So, the GALI method allows us to study more efficiently 
the geometrical properties of the dynamics in the neighbor- 
hood of an orbit, especially in higher dimensions, where it 
allows for a much faster determination of its chaotic nature, 
overcoming the limitations of the SALI method. In the case 
of regular motion, GALIk is either a constant, or decays by 
power laws that depend on the dimensionality of the sub- 
space in which the orbit lies, which can prove useful e.g., if 
our orbits are in a "sticky" region, or if our system happens 
to possess fewer or more than independent integrals of the 
motion (i.e. is partially integrable or super-integrable respec- 
tively). It should be mentioned, that both the SALI and the 
GALI methods are very reliable and efficient. However, in 
order to apply them we have to use the variational equations 
and also several sets of deviation vectors. On the contrary, 
the dFNVI method is much faster than these methods be- 
cause it does nor require the use of the variational equations 
and the deviation vectors. Simply we follow the time evolu- 
tion of an orbit with initial condition x(0) and computing its 
norm by integrating only the basic equations of motion. 

In what follows, we shall present and compare repre- 
sentative examples for 2D and 3D orbits of the dynamical 
systems (7) and (6) respectively, using the above four indi- 
cators. In Figure 16a-d, we can see the results provided by 
the four indicators for a regular 2D orbit with initial con- 
ditions: xo = 0.27, Jo = and = 0, while the initial 
value of /7vo is found from the energy integral (7). The value 
of the energy is hi - 1. In Fig. 16a we present the evolu- 
tion of the DFLI for a time period of 5 x 10^ time units. We 



observe that the DFLI fluctuates around small values (DFLI 
< 10), which indicates regular motion. The initial deviation 
vector, w = {dx, dy, dp^, dpy), used for the orbit shown in 
Fig. 16a is w(0) - (1, 0, 0, 0). Fig. 16b depicts the evolution 
of the RLI for a time interval of 10^ time units, for the same 
regular orbit. The RLI clearly reveals the ordered nature of 
the orbit, since RLI < 10"^'^ throughout the time evolution. 
The initial deviation vectors, used for the computation of the 
RLI are wi(0) = (1,0,0,0) and W2(0) = (1,0,0,0), while 
6x = 10"'^. Note that the t axis is in log scale. In Fig. 16c we 
present the evolution of SALI for a time interval of 5 x 10^ 
time units. Also this method indicates the regular character 
of the orbit, since it exhibits small fluctuations around non- 
zero values. The initial deviation vectors used for the orbit 
of Fig. 16c are wi(0) = (1,0,0,0) and W2(0) = (0,0, 1,0), 
but in general, any other choice of initial values leads to sim- 
ilar behavior of the SALI. From Eqs. (18) and (19) it follows 
that in the case of a Hamiltonian system with = 2 degrees 
of freedom GALI2 wiU always remains different from zero, 
while GALI3 and GALI4 should decay to zero following a 
power law, whose exponent depends on the number m of 
the deviation vectors that are initially tangent to the torus 
(19^ which the orbit lies. Now, for the regular orbit of the 2D 
Hamiltonian (7) and for a random choice of initial deviation 
vectors, we expect the GALIk indices to behave as 

GALI2 (f) oc constant, 
GALI3 (t) oc 1 



GALI4 (0 oc 



1 



(20) 



Fig. 16d shows the behavior of GALl^ for the same regular 
2D orbit. Indeed, the GALlk indices obey to the power laws 
respectively given in Eq. (20). In Fig. 16d with red color we 
plot the power law 1 /f^, while the power law 1 /f^ is colored 
in green. Note that the t axis is in log scale. The evolution of 
the FN VI for this regular 2D orbit, for a time interval of 10^ 
time units, is presented in Fig. 3a, while dFNVI = 0.009. 
We see, that all the outcomes derived using five different 
dynamical indicators coincide to the ordered nature of the 
orbit. 

Figure 17a-d is similar to Fig. 16a-d but for a chaotic 
2D orbit. In Fig. 17a-d, we present the results provided by 
the four dynamical indicators for a chaotic 2D orbit with 
initial conditions: xq = 0.82, >'o = and pxo = 0, while 
the initial value of Pyo is found from the energy integral 
(7). The value of the energy is hi = 1. In Fig. 17a we 
present the evolution of the DFLI for a time period of 5 X 
lO"' time units. We observe that the DFLI displays an expo- 
nential growth with time, which indicates chaotic motion. 
The initial deviation vector used for the orbit shown in Fig. 
17a is w(0) = (1,0,0,0). Fig. 17b depicts the evolution of 
the RLI for a time interval of 10^ time units, for the same 
chaotic orbit. Approximately for the first 350 time units, the 
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RLI has low values (RLI < 10"^) but for the rest time inter- 
val RLI > 10 Note that the t axis is in log scale. There- 
fore, we may conclude that the orbit is chaotic. Note, that 
the RLI method gives inconclusive or even misleading re- 
sults for very small time intervals of the numerical integra- 
tion. The initial deviation vectors, used for the computation 
of the RLI are Wi(0) = (1,0,0,0) and W2(0) = (1,0,0,0), 
while 6x - 10"'^. In Fig. 17c we present the evolution of 
the SALI. We see that after a small transient time, the SALI 
falls abruptly to zero. At f ^ 450 time units the SALI be- 
comes zero, as it has reached the limit of the accuracy of 
the computer (10""'), which means that the two deviation 
vectors have the same direction. Thus, after f ^ 450 time 
units the two normalized vectors are represented by exactly 
the same numbers in the computer and we can safely argue, 
that to this accuracy the orbit is chaotic. The initial deviation 
vectors used for the orbit of Fig. 17c are Wi(0) - (1,0,0,0) 
and W2(0) - (0,0, 1,0). From Eq. (17) it is evident that in 
the case of a Hamiltonian system with N - 2 degrees of 
freedom, GALI^ indices should decrease exponentially ap- 
proaching the zero-value. In Fig. 17d we observe the be- 
havior of GALIk for the same chaotic 2D orbit. Indeed, the 
GALIk indices obey to the exponential law given in Eq. (17). 
The evolution of the FNVI for this chaotic 2D orbit, for a 



time interval of 10^ time units, is presented in Fig. 3b, while 
dFNVI - 0.45. We see, that all the outcomes obtained using 
five different dynamical indicators coincide to the chaotic 
nature of the orbit. 

We shall now proceed, presenting and comparing two 
more examples regarding 3D orbits of the Hamiltonian sys- 
tem (6). In Figure 18a-d, we can see the results provided by 
the four indicators for a regular 3D orbit with initial condi- 
tions: X() - 0.09, yo = 0,zo = 0.1,pj.o - PzO - 0, while the 
initial value of pyQ is found from the energy integral (6). The 
value of the energy is hj, - 1 . In Fig. 18a we present the evo- 
lution of the DELI for a time period of 5 x lO'' time units. We 
observe that the DELI fluctuates around small values (DELI 
< 10), which indicates regular motion. The initial deviation 
vector, w = {dx,dy,dz,dpx,dpy,dpS), used for the orbit 
shown in Fig. 18a is w(0) = (1,0,0,0,0,0). Fig. 18b de- 
picts the evolution of the RLI for a time interval of 10^ time 
units, for the same regular orbit. The RLI clearly reveals the 
ordered character of the orbit, since RLI < 10~^2 through- 
out the time evolution. Note that the t axis is in log scale. 
The initial deviation vectors, used for the computation of the 
RLI are wi(0) = (1, 0, 0, 0, 0, 0) and W2(0) = (1, 0, 0, 0, 0, 0), 
while 6x - 10"^^^. In Fig. 18c we present the evolution of 
SALI for a time interval of 5 x 10^ time units. Also this 
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Fig. 17 (a-d): Similar to Fig. 16 a-d but for a chaotic 2D orbit of the dynamical system (7). The initial conditions of the 2D orbit and more details 
are provided in the text. 



method indicates the regular character of the orbit, since it 
exhibits small fluctuations around non-zero values (SALI 
e [1, V2]). The initial deviation vectors used for the or- 
bit of Fig. 18c are wi(0) = (1,0,0,0,0,0) and W2(0) = 
(0, 0, 0, 1, 0, 0), but in general any other choice of initial val- 
ues leads to similar behavior of the SALI. From Eqs. (18) 
and (19) it follows that in the case of a Hamiltonian system 
with A'^ = 3 degrees of freedom GALI2 and GALI3 will al- 
ways remain different from zero, while GALI4, GALI5 and 
GALIfi should decay to zero following a power law, whose 
exponent depends on the number m of deviation vectors that 
are initially tangent to the torus on which the orbit lies. Now, 
for the regular orbit of the 3D Hamiltonian (6) and for a ran- 
dom choice of initial deviation vectors, we expect the GALlk 
indices to behave as 
GALI2 (t) oc constant, 
GALI3 (t) oc constant, 

GALI4 it) oc 



1 



GALI5 it)c^^, 
GALlg (0 « ^■ 



(21) 



Fig. 18d shows the behavior of GALlk for the same regular 
3D orbit. Indeed, the GALlk indices obey to the power laws 



respectively given in Eq. (21). In Fig. 18d with red color we 
plot the power law l/t^, the power law 1/?'* is colored in 
green, while the power law 1 /t^ is plotted with blue color. 
Note that the t axis is in log scale. With a more closer look in 
Fig. 18d, we observe that there is a superposition of GALI2 
and GALI3 as they evolve almost identically. The evolution 
of the FNVI for this regular 3D orbit, for a time interval of 
10^ time units, is presented in Fig. 4a, while dFNVI = 0.006. 
We see, that once more, all the outcomes derived using five 
different dynamical indicators coincide to the ordered nature 
of the orbit. 

Finally, Figure 19a-d is similar to Fig. 18a-d but for a 
chaotic 3D orbit. In Fig. 19a-d, we present the results pro- 
vided by the four dynamical indicators for a chaotic 3D or- 
bit with initial conditions: xq = 0.86, = 0,zo = 0.1,^^0 = 
p^o - 0, while the initial value of pyo is found from the en- 
ergy integral (6). The value of the energy is hs = 1. In Fig. 
19a we present the evolution of the DFLI for a time period of 
5 X lO-' time units. We observe that the DFLI displays a very 
rapid growth with time, which indicates chaotic motion. The 
initial deviation vector used for the orbit shown in Fig. 19a is 
w(0) = (1,0, 0, 0). Fig. 19b depicts the evolution of the RLI 
for a time interval of 10^ time units, for the same chaotic or- 
bit. Approximately for the first 200 time units, the RLI has 
low values (RLI < 10"^) but for the rest time interval RLI 
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> 10" . Note that the f axis is in log scale. Therefore, we 
may conclude that the orbit is chaotic. Note, that again the 
RLl method gives inconclusive or even misleading results 
for very small time intervals of the numerical integration. 
The initial deviation vectors, used for the computation of the 
RLI are Wi(0) = (1,0, 0, 0, 0, 0) and W2(0) = (1,0, 0, 0, 0, 0), 
while Sx - 10 In Fig. 19c we present the evolution of 
the SALl. We see that the SALl decreases very rapidly ap- 
proaching the zero-value. At f a; 250 time units the SALl be- 
comes zero, as it has reached the limit of the accuracy of the 
computer (10""'), which means that the two deviation vec- 
tors have the same direction. Thus, after t ^ 250 time units 
the two normalized vectors are represented by exactly the 
same numbers in the computer and we can safely argue, that 
to this accuracy the orbit is chaotic. The initial deviation vec- 
tors used for the orbit of Fig. 19c are Wi(0) = (1,0, 0, 0, 0, 0) 
and W2(0) = (0, 0, 0, 1, 0, 0). From Eq. (17) it is evident that 
in the case of a Hamiltonian system with N = 3 degrees of 
freedom, all GALl^ indices should decrease exponentially 
approaching the zero-value. In Fig. 19d we observe the be- 
havior of GALIk for the same chaotic 3D orbit. Indeed, all 
the GALIk indices obey to the theoretically predicted law 
given in Eq. (17). The evolution of the FNVI for this chaotic 
3D orbit, for a time interval of 10^ time units, is presented in 



Fig. 4b, while dFNVI - 0.52. We see, that all the outcomes 
obtained using five different dynamical indicators coincide 
to the chaotic nature of the orbit. 

From the above examples, it becomes evident that all the 
different methods we used (DFLI, RLI, SALl and GALI) 
need not only the computation of the basic equations of mo- 
tion but also the simultaneous computation of the so-called 
variational equations and of course several sets of deviation 
vectors. As the total number of the equations and inevitably 
the total number of the variables increases, we have also an 
increase to the time needed from the computational code in 
order to integrate the system of the differential equations and 
provide the results. Thus, our proposed FNVI method has 
one important advantage over all the other methods men- 
tioned above, since it requires only the computation of the 
basic set of the equations of motion without any use of vari- 
ational equations and deviation vectors. Let us present a spe- 
cific illuminating example regarding the speed of the meth- 
ods. Figure 20 presents the CPU time in h for each dynam- 
ical method required, using the same technique described 
in section 3 and of course the same processor, in order to 
construct the grid-plot shown in Fig. 7a. We observe, that 
the dFNVI method, needs only about 7 h of CPU time. The 
total time required by the S(g) spectrum is similar to the 
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Fig. 19 (a-d): Similar to Fig. 18 a-d but for a chaotic 3D orbit of the dynamical system (6). The initial conditions of the 3D orbit and more details 
are provided in the text. 



time needed for the dFN VI method. On the other hand, if we 
choose to use methods such as the DFLI, SALI and GALI^, 
that need the additional computation of the variational equa- 
tions, the total time required for the completion of the com- 
putational task increases significantly. Furthermore, we see 
from Fig. 20 that the most time consuming methods are the 
LCE and RLI. This is true, because when we use these two 
methods, we have to integrate each orbit at least for 10'' time 
units in order to obtain reliable information about the or- 
dered or chaotic nature of the orbit. Thus, it becomes more 
than obvious, than the use of these two methods doubles the 
total CPU time needed by the dFNVI method. Therefore, 
our new method is not only reliable but also very fast com- 
pared with other methods and can be applied easily when we 
need to integrate a large number of orbits so as to construct 
a grid-plot of initial conditions (see Figs. 7 and 10). 

Apart from the dynamical indicators described and ap- 
phed previously, there are also other methods for distinguish- 
ing between ordered and chaotic motion, such as the mean 
exponential growth of nearby orbits (MEGNO) [4, 5], the 
power spectra of deviation vectors [30], the method of the 
low frequency power (LFP) [13, 29], the "0-1" test [10], as 
well as some other more recently introduced techniques [11, 
24]. Moreover, a method that have also the ability to iden- 
tify sticky orbits is the FtDt method [12]. This method uses 
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Fig. 20 The total CPU time in h for each dynamical method required, 
in order to construct the grid-plot shown in Fig. 7a. See text for more 
details. 



the Fast Fourier Transform of a series of time interval, each 
one representing the time that elapsed between two succes- 
sive points on the Poincare surface of section. However, we 
do not feel that it is necessary to provide examples using all 
of these methods. We beUeve, that the comparison between 
the FNVI method and four well-known dynamical indica- 
tors proved its reliabihty and also revealed its remarkable 
efficiency. 



Detecting order and chaos in Hamiltonian systems by the Fast Norm Vector Indicator (FNVI) method 



23 



6 Discussion and conclusions 

In the present paper, we have tried to introduce a new, fast, 
efficient and easy to compute method in order to distinguish 
between order and chaos in 2D and 3D autonomous Hamil- 
tonian systems. We have also conducted a detailed study of 
the behavior of this new indicator for both chaotic and reg- 
ular orbits. The main results of this research can be summa- 
rized as follows: 

1. The FNVI method proves to be an ideal detector of 
chaoticity independent of the dimensions of the dynamical 
system. It displays large, abrupt and random fluctuations for 
chaotic orbits, while it remains almost constant for ordered 
ones and so it clearly distinguishes between these two cases. 
Its main advantages are its simplicity, efficiency and relia- 
bility as it can rapidly and accurately determine the chaotic 
versus ordered nature of a given orbit. Of course, the FNVI 
method can provide only qualitative results regarding the na- 
ture of an orbit and therefore, we have to inspect the shape of 
FNVI each time in order to characterize an orbit. Obviously, 
this is not very practical when someone wants to check a 
large volume of orbits, so as to form an idea about the global 
structure of the dynamical system. Therefore, we proceeded 
one step further estabhshing a numerical criterion in order 
to quantify ffie results obtained by the FNVI method. This 
criterion derived by exploiting the shape of the FNVI. When 
the orbit is regular the FNVI remains almost constant, while 
in the case of a chaotic orbit it displays high fluctuations. 
Thus, we calculate the maximum and the minimum value of 
FNVI when t e [200, 1000]. We choose this particular time 
interval because in the case where the orbit is regular, for 
t < 200 time units there is a transient period of fluctuation, 
which may cause a problem to our criterion. Using the above 
procedure we defined the dFNVI. The value of the dFNVI 
can provide us the quantitative criterion that we seek. We 
point out, that it is not very easy to define a threshold value, 
so ffiat the dFNVI being larger ffian this value reliably sig- 
nifies chaoticity. Nevertheless, extensive numerical experi- 
ments in both dynamical systems indicate that in general, a 
good guess for this value could be 0.05. Thus, when dFNVI 
> 0.05 the orbit is chaotic, while when dFNVI < 0.05 the 
orbit ir ordered. This threshold value applies in both 2D and 
3D dynamical systems. 

2. We emphasize, that the main advantage of the FNVI 
method, is that it needs only the computation of the basic 
equations of motion. In particular, we only have to follow 
the evolution of the norm of the vector x(f) and compare 
its value in every time step throughout the entire time inter- 
val of the numerical integration with the initial value x(0) 
according to Eq. (4). The FNVI method is very fast, as it 
requires only lO"' time units of numerical integration, in or- 
der to provide reliable and conclusive results regarding the 
character of an orbit. Of course, there are also other meth- 



ods and indicators capable to discriminate between ordered 
and chaotic orbits. However, the vast majority of them need 
the computation of the variational equations and inevitably 
the use of several sets of initial deviation vectors in order 
to function and provide their results. The additional use and 
therefore the computation of more equations increases sig- 
nificantly the total time needed from the computational rou- 
tine to provide the output in each case (see Fig. 20). Thus, 
our proposed FNVI method has an important advantage over 
most the other dynamical methods, since it requires only 
the computation of the basic set of the equations of motion. 
Therefore, our method is very fast which makes it is an ideal 
tool when we need to calculate a large number of orbits so 
as to construct a grid of initial conditions. 

3. Using the numerical quantified version of FNVI, that 

is the dFNVI, we are in a position to characterize reliable an 
orbit of being chaotic or ordered. Exploiting the advantages 
of the dFNVI method, we have constructed detailed phase- 
space portraits (grid-plots) both for 2D and 3D Hamiltonian 
systems, where the chaotic and ordered regions are clearly 
distinguished (see Figs. 7 and 10). We were thus able to 
trace in a fast and systematic way very small islands of or- 
dered motion, whose detection by traditional methods would 
be very difficult, if not impossible, and moreover very time 
consuming. This approach is therefore, expected to provide 
useful tools for the location of stable periodic orbits, or the 
computation of the phase-space volume occupied by ordered 
or chaotic motion in multi-dimensional systems (N > 3), 
where the PSS plots can not be easily visualized and inter- 
preted and furthermore, very few other similar techniques 
of practical value are available. Here, we have to point out 
that grid-plots can be also constructed using other dynam- 
ical indicators. For the reasons mentioned in the previous 
point, the dFNVI method is much more faster than most of 
the other methods. 

4. We used the new FNVI method in order to follow 

the time evolution of a sticky orbit. The study of the evo- 
lution of sticky orbits (orbits that, although they are chaotic, 
they are restricted in a thin layer of phase space for a large 
value of integration time) is crucial for the understanding 
of the structure of the dynamical system. The ability of our 
new method to extract reliable information concerning the 
chaotic or regular behavior of an orbit enables us to over- 
come the ambiguity in the case of sticky orbits. In the case 
of sticky orbits, the distinction can not be easily made using 
the PSS technique, as a sticky orbit can be mistaken for an 
ordered one, if not enough time is given and the orbit has 
not yet escaped to the surrounding chaotic domain (see Fig. 
13a-d). Moreover, the LCE is not sensitive enough to iden- 
tify the difference between a regular and a sticky orbit (see 
Fig. 15a). On the contrary, the FNVI method can be applied 
to obtain a clear and reliable distinction between ordered and 
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Sticky orbits (see Fig. 15b) and this could be regarded as a 
major advantage of our method. 

In the present research, we have tested approximately 
10^ orbits (2D and 3D) in the dynamical systems (6) and 
(7). Since our results are consistent, we may safely con- 
clude that the FNVI method (and also its numerical version 
the dFNVI) is a very fast and reliable tool for distinguish- 
ing between order and chaos in Hamiltonian systems. It is 
om future plans, to apply this new method in other kinds of 
interesting potentials and also in more complicated dynam- 
ical systems, in order to test further its efficiency and reli- 
abihty. Furthermore, we will apply this new method in or- 
der to study the evolution of three-dimensional sticky orbits. 
Moreover, with additional theoretical work, we will inves- 
tigate if it always discriminates correctly between ordered 
and chaotic motion in Hamiltonian systems, especially in 
dynamical systems with three degrees of freedom. 
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